# Analyzing the results of LS-pfOCE results over full 3D OCE measurmeent 
# output scatter plot indicating stiffness gradients of LS-pfOCE 

# Created by Yuechuan Lin 
# 2021-04-13 
from numpy.core.arrayprint import _void_scalar_repr
from numpy.lib.recfunctions import _drop_fields_dispatcher
import matplotlib.mlab as mlab
import pandas as pd 
import os
from PyOCT import misc 
import numpy as np 
import matplotlib.pyplot as plt 
import matplotlib 
from matplotlib import cm 
import matplotlib.ticker as tck
import matplotlib.ticker 
import h5py 
import time
from scipy import ndimage
from scipy.io import loadmat
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,
                               AutoMinorLocator)
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.collections import PathCollection
from sklearn.neighbors import KernelDensity
import statsmodels.api as sm
# set font of plot 
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = ['Helvetica']
font = {'weight': 'normal',
        'size'   : 10}
matplotlib.rc('font', **font)
matplotlib.rc('text', usetex=False)
matplotlib.rcParams['pdf.fonttype'] = 42
matplotlib.rcParams['ps.fonttype'] = 42
matplotlib.rcParams['mathtext.default'] = 'regular'

def makeColors(vals,cmap="plasma",alpha=0.5):
    colors = np.zeros((len(vals),3))
    print("Colormap kde with min {} and max {}".format(vals.min(),0.7*vals.max()))
    norm = matplotlib.colors.Normalize(vmin=vals.min(),vmax=0.7*vals.max())
    colors = [cm.ScalarMappable(norm=norm,cmap=cmap).to_rgba(val,alpha=alpha) for val in vals]

    vals = vals*1e6
    cml = cm.ScalarMappable(norm=matplotlib.colors.Normalize(vmin=vals.min(),vmax=0.7*vals.max()),cmap=cmap)
    max_v = 0.7*vals.max()
    return colors, cml, max_v 

def fmt(x, pos):
    #a, b = '{:.1e}'.format(x).split('e')
    #b = int(b)
   # x = x*1e6
    a, b = '{:.1f}'.format(x).split('.')
    #return r'${} \times 10^{{{}}}$'.format(a, b)
    return a

def LoadBeadWiseData(fileName,xoffset,cellNum):
    """
    Load measured baed-wise results, which are in a format of:
    [z x y snr ampRatio stdA stdphi Atot phitot APT phiPT Amech phimech Frad GRe GIm]
    Input:
    : fileName: including the data path and file suffix. 
    : xoffset: the offset between the cetners of lightsheet force to OCT FOV 
    """
    fid_beads = h5py.File(fileName,"r")
    results = np.transpose(fid_beads["results_beadwise"][()]) #(N,16) 
    mask = np.any(np.isnan(results),axis=1) 
    results = results[~mask] 
    fid_beads.close()

    #now exclue the depth limit outside [-40,40]um along z 
    results_tmp = [] 
    for i in range(np.shape(results)[0]):
        if np.abs(results[i,0]) < 40: 
            if (results[i,14] > 30) and (results[i,15]>0):
                results_tmp.append(results[i,:]) 
    results_output = np.asarray(results_tmp) 

    z = results_output[:,0]
    x = results_output[:,1] - xoffset
    y = results_output[:,2] 
    Gre = results_output[:,14]
    Gim = results_output[:,15] 
    Gmag = np.sqrt(Gre**2+Gim**2) 
    R = Gim/Gre 
    snr = results_output[:,3] 
    r = np.zeros(np.shape(z))
    group = np.zeros(np.shape(z),dtype=int)
    group_rmedian = np.zeros(np.shape(z)) 
    cell_num_mat = cellNum*np.ones(np.shape(z))
    df_res = pd.DataFrame({"cellNum":cell_num_mat,"z":z,"x":x,"y":y,"Gre":Gre,"Gim":Gim,"Gmag":Gmag,"R":R,"snr":snr,"r":r,"group":group,"group_rmedian":group_rmedian})
    return df_res 

def LoadCellVol(fileName,xoffset):
    fid_cell = loadmat(fileName) 
    cellbody = fid_cell["cell"] 
    zc = cellbody[:,0]
    xc = cellbody[:,1] - xoffset
    yc = cellbody[:,2] 
    #stdc = cellbody[:,3]
    intc = cellbody[:,4] 
    return zc,xc,yc,intc

def LoadingFromFileList(root_path,beadsFileNames,cellFileNames,xoffset,CenterToCell=True,normBkgnd=False):
    """
    Loading all beads and cell data file. 
    Input:
    : root_path: path where all data files saved 
    : beadsFileNames, cellFileNames: list, data path to beads, cell data respectively. 
    : xoffset: input Ndarray, indicating the xoffset for each cell data 
    : CenterToCell: making the mass center of cell as the coordinate center of all data 
    : normBkgnd: normalizing the output results of mechancis to the background value. 
    """
    N_cell = len(beadsFileNames) 
    df_res = pd.DataFrame(columns=["cellNum","z","x","y","Gre","Gim","Gmag","R","snr","r","group","group_rmedian"])
    for i in range(N_cell):
        print("Loading {} of {} cells...".format(i+1,N_cell))
        if beadsFileNames[i].endswith(".mat"):
            tmpBeadsFileName = beadsFileNames[i] 
        else:
            tmpBeadsFileName = beadsFileNames[i] + ".mat" 
        if cellFileNames[i].endswith(".mat"):
            tmpCellFileName = cellFileNames[i] 
        else:
            tmpCellFileName = cellFileNames[i] + ".mat"
        tmp_df_res = LoadBeadWiseData(os.path.join(root_path,tmpBeadsFileName),xoffset[i],i) 
        tmpzc,tmpxc,tmpyc,tmpIntc = LoadCellVol(os.path.join(root_path,tmpCellFileName),xoffset[i])
        [zcen,xcen,ycen] = [np.sum(tmpzc*tmpIntc)/np.sum(tmpIntc),np.sum(tmpxc*tmpIntc)/np.sum(tmpIntc),np.sum(tmpyc*tmpIntc)/np.sum(tmpIntc)]
        if CenterToCell:
            tmp_df_res["z"] = tmp_df_res["z"] - zcen 
            tmp_df_res["x"] = tmp_df_res["x"] - xcen 
            tmp_df_res["y"] = tmp_df_res["y"] - ycen 
            tmpzc = tmpzc - zcen 
            tmpyc = tmpyc - ycen 
            tmpxc = tmpxc - xcen 
        for k in range(len(tmp_df_res.index)):
            tmprArr = np.sqrt((tmp_df_res.loc[k,"z"]-tmpzc)**2+(tmp_df_res.loc[k,"x"]-tmpxc)**2+(tmp_df_res.loc[k,"y"]-tmpyc)**2)
            tmprArr = np.sort(tmprArr) 
            tmp_df_res.loc[k,"r"] = np.amin(tmprArr) #tmprArr[np.argmin(tmprArr)]
        if normBkgnd:
            tmp_bkgnd_df = tmp_df_res[tmp_df_res["r"]>np.percentile(tmp_df_res["r"],90)] 
            print("Num of bkgnd indx is {}".format(len(tmp_bkgnd_df.index))) 
            print("Gre bkgnd is : {} Pa".format(np.median(tmp_bkgnd_df["Gre"])))
            print("Gim bkgnd is : {} Pa".format(np.median(tmp_bkgnd_df["Gim"])))
            print("R bkgnd is {}".format(np.median(tmp_bkgnd_df["R"])))
            tmp_df_res["Gre"] = tmp_df_res["Gre"]/np.median(tmp_bkgnd_df["Gre"])
            tmp_df_res["Gim"] = tmp_df_res["Gim"]/np.median(tmp_bkgnd_df["Gim"])
            tmp_df_res["Gmag"] = tmp_df_res["Gmag"]/np.median(tmp_bkgnd_df["Gmag"])
            tmp_df_res["R"] = tmp_df_res["R"]/np.median(tmp_bkgnd_df["R"])
        df_res = df_res.append(tmp_df_res,ignore_index=True)   
    return df_res 

def SelectROI(ROIList,df_input,perRange=[10,90],normBkgnd=False):
    """
    Group beads-wise data results based on the distance w.r.t cell body. 
    Input:
    : ROIList: list of ROI as dimension of [N,2] where N indicates the number of ROI, 
        first column is the minimum distance in the group while second column is the largest distance. 
    : df_input: beads-wise DataFrame results 
    : perRange: percentile, only the results in the range of percentile larger than perRange[0] and smaller than perRange[1] will be included. 
    : normBkgnd: normalizing all mechanics to bkgnd. 
    """
    ROIList = np.asarray(ROIList)
    df_input["group"] = -1 
    df_input["group_rmedian"] = -1 
    df_input["rGmean"] = -1 
    df_input["rGcon"] = -1 
    df_input["rRmean"] = -1 
    df_input["rRcon"] = -1 
    df_input["rrmean"] = -1 
    df_input["rrcon"] = -1 
    #assign group 
    for i in range(np.shape(ROIList)[0]):
        df_input.loc[df_input.index[(df_input["r"]>ROIList[i,0]) & (df_input["r"]<=ROIList[i,1])],"group"] = int(i) 
   # df_input.loc[df_input.index[(df_input["r"] > ROIList[-1,1])],"group"] = int(np.shape(ROIList)[0]) 
    
    #only accepth the values in the range of perRange within each group
    for i in range(np.shape(ROIList)[0]):
        tmp_min = np.percentile(df_input.loc[df_input["group"] == i,"Gre"],perRange[0]) 
        tmp_max = np.percentile(df_input.loc[df_input["group"] == i,"Gre"],perRange[1]) 
        df_input.loc[(df_input["group"] == i) & (df_input["Gre"]>tmp_max),"group"] = -1 
        df_input.loc[(df_input["group"] == i) & (df_input["Gre"]<tmp_min),"group"] = -1 
        
    #calculate the position median 
    for i in range(np.shape(ROIList)[0]):
        print("ROI [{},{}] has {} beads".format(ROIList[i,0],ROIList[i,1],len(df_input[df_input["group"] == i].index)))
        df_input.loc[df_input["group"] == i,"group_rmedian"] =  np.median(df_input.loc[df_input["group"] == i,"r"])
        tmpN = len(df_input.index[df_input["group"] == i])
        tmp_rGmean, tmp_rGcon = misc.mean_confidence_interval(df_input.loc[df_input["group"] == i,"Gre"],confidence=0.95,inputType="d")
        tmp_rRmean, tmp_rRcon = misc.mean_confidence_interval(df_input.loc[df_input["group"] == i,"R"],confidence=0.95,inputType="d") 
        tmp_rrmean, tmp_rrcon = misc.mean_confidence_interval(df_input.loc[df_input["group"] == i,"r"],confidence=0.95,inputType="d") 
        df_input.loc[df_input["group"] == i,"rGmean"] = np.ones((tmpN,)) * tmp_rGmean
        df_input.loc[df_input["group"] == i,"rGcon"] = np.ones((tmpN,)) * tmp_rGcon
        df_input.loc[df_input["group"] == i,"rRmean"] = np.ones((tmpN,)) * tmp_rRmean
        df_input.loc[df_input["group"] == i,"rRcon"] = np.ones((tmpN,)) * tmp_rRcon
        df_input.loc[df_input["group"] == i,"rrmean"] = np.ones((tmpN,)) * tmp_rrmean
        df_input.loc[df_input["group"] == i,"rrcon"] = np.ones((tmpN,)) * tmp_rrcon

    return df_input 

if __name__ == "__main__":
    colors = ["#019BD8", "#D81C28","#FF7F0E","#2CA02C","#e377c2"]
    scatterSize = 15
    lw = 1.6
    label_fontsize=12 
    logScale = True 
    plotR = True 
    normBkgnd = True 
    if normBkgnd:
        curveFit = True 
    else:
        curveFit = True 
    root_path =  "./Fig3_data"   
    save_path = "" #define your savepath here
    if not os.path.isdir(save_path):
        os.mkdir(save_path) 
    
    NormalOCE_beadsFileNames = ["earlyCell11/eCell11/BeadsRes_210128","earlyCell3/eCell3/BeadsRes_201224","earlyCell4/eCell4/BeadsRes_210107","earlyCell6/eCell6/BeadsRes_210110","earlyCell7/eCell7/BeadsRes_210111","earlyCell8/eCell8/BeadsRes_210123","earlyCell9/eCell9/BeadsRes_210124","earlyCell10\eCell10/BeadsRes_210126"] #'/CellI/beadswise_CellI_200821.mat','/CellII/beadswise_CellII_200821.mat','/CellIII/beadswise_CellIII_200911.mat','/CellIV/beadswise_CellIV_200923.mat',"/CellV/beforeCytoD/beadswise_CellV_201012.mat",'/CellVI/beforeCytoD/beadswise_CellVI_201027.mat','/CellVII/beforeCytoD/beadswise_CellVII_201027.mat','/postCellI/beforeCytoD/beadswise_postCellI_201123.mat']
    NormalOCE_cellFileNames = ["earlyCell11/eCell11/CellBody_210128","earlyCell3/eCell3/CellBody_201224","earlyCell4/eCell4/CellBody_210107","earlyCell6/eCell6/CellBody_210110","earlyCell7/eCell7/CellBody_210111","earlyCell8/eCell8/CellBody_210123","earlyCell9/eCell9/CellBody_210124","earlyCell10\eCell10/CellBody_210126"] #'/CellI/CellI_200821.mat','/CellII/CellII_200821.mat','/CellIII/CellIII_200911.mat','/CellIV/CellIV_200923.mat','/CellV/beforeCytoD/CellV_201012.mat','/CellVI/beforeCytoD/CellVI_201027.mat','/CellVII/beforeCytoD/CellVII_201027.mat','/postCellI/beforeCytoD/postCellI_201123.mat']
    
    NormalOCE_N_cell = len(NormalOCE_beadsFileNames)
    NormalOCE_xoffset = np.asarray([52,52,52,52,52,52,52,52,52,52,52,52,52,52,52]) 

    NormalOCE_ROI = np.asarray([[0.0,3.0],[6.1,10],[11.1,14],[18.1,21],[25,27],[30,35],[40,45],[50,55],[60,65],[70,85],[90,105],[106,115]])
    df_res = LoadingFromFileList(root_path,NormalOCE_beadsFileNames,NormalOCE_cellFileNames,NormalOCE_xoffset,CenterToCell=True,normBkgnd=normBkgnd) 
    df_res = SelectROI(NormalOCE_ROI,df_res,perRange=[10,95])
    df_res.to_csv(os.path.join(save_path,"results.csv"),index=False) 
    print("Normal Cell OCE results...")
    print(df_res)

    #plot Gre w.r.t dist_nearest 
    figK,axK = misc.create_fig(figsize=[4,2],bottom=0.2,left=0.15,right=0.9,top=0.8)
    axK.set_xlabel(r"Distance to cell $r$ ($\mathrm{\mu}$m)",fontsize=label_fontsize,labelpad=0)
    axK.set_ylabel(r"Elasticity $\mathrm{G^{\prime}}$ (Pa)",fontsize=label_fontsize,labelpad=0) 
    Kde_n = KernelDensity(kernel='exponential', bandwidth=5,metric="chebyshev")
    r_norm = df_res["r"].to_numpy()
    Gre_norm = df_res["Gre"].to_numpy() 
    r_norm = np.reshape(r_norm,(np.size(r_norm),1)) 
    Gre_norm = np.reshape(Gre_norm,(np.size(Gre_norm),1)) 
    normX =  np.concatenate((r_norm,Gre_norm),axis=1) 
    Kde_n.fit(normX)  
    distributed_normGre = np.exp(Kde_n.score_samples(normX)) 
    normColors,ncml,nmax = makeColors(distributed_normGre,cmap="viridis",alpha=0.7)
    axK.scatter(df_res["r"],df_res["Gre"],color=normColors,s=10,rasterized=True) 
    ncbar = plt.colorbar(ncml,ax=axK,fraction=0.10,shrink=1.0,aspect=10,format=matplotlib.ticker.FuncFormatter(fmt),ticks=[25,50,75,100,nmax])
    ncbar.ax.set_title(r"p.d.($\times 10^{-6}$)",fontsize=label_fontsize)
    ncbar.ax.set_yticklabels(['25', '50', '75','100','>'+str(int(nmax))])
    axK.set_ylim([0,400])
    axK.set_yticks([0,100,200,300,400])   


    #plot the trend 
    fig1, ax1 = misc.create_fig(figsize=[3.2,2.7],top=0.9,bottom=0.2,left=0.2,right=0.84) 
    if curveFit:
        axCF = fig1.add_axes([0.45,0.6,0.38,0.292])
        axCF.set_xscale("log")
        axCF.set_yscale("log")
        axCF.tick_params(axis='both', labelsize=5,pad=0)
        axCF.set_yticks([], minor=True)
        axCF.set_yticks([1.0,1.5,2.0,2.5])
        axCF.set_yticklabels(["1.0","1.5","2.0","2.5"],fontsize=7)
        axCF.set_xlim([0.6,60])
        axCF.set_ylim([0.9,3.0])
        axCF.set_xlabel(r"r",fontsize=7,labelpad=0.5)
        axCF.set_ylabel(r"$G^{\prime}$",fontsize=7,labelpad=0.5)


    if normBkgnd: 
        ax1.set_ylabel(r"Rel. Elasticity $\mathrm{G^{\prime}}$",labelpad=0,fontsize=label_fontsize) 
    else:
        if logScale:
            ax1.set_ylabel(r"Elasticity $\mathrm{G^{\prime}}$ ($\times 10^2$ Pa)",labelpad=0,fontsize=label_fontsize)
            ax1.set_yscale("log")
        else:
            ax1.set_ylabel(r"Elasticity $\mathrm{G^{\prime}} (Pa)$",labelpad=0,fontsize=label_fontsize) 
        
    if plotR:
        ax11 = ax1.twinx()
        #if logScale:
        #    ax11.set_yscale("log")
        ax1.set_zorder(ax11.get_zorder()+1) 
        ax1.patch.set_visible(False) 
        ax11.tick_params(axis="y",which="both",left=False,labelleft=False,length=2.0,right=True,labelright=True) 
        ax11.set_ylabel(r"Rel. viscosity $\mathrm{R}$", labelpad=0,fontsize=label_fontsize) 
        if not normBkgnd:
            ax11.set_ylim([0.01,5.0])
            ax11.set_yticks([0.5,1.0,2.0,3.0,4.0,5.0]) #ax11.set_yticks([0.1,1,2,3,4])
            ax11.set_yticklabels([])
            ax11.set_yticklabels(["0.5","1.0","2.0","3.0","4.0","5.0"])
           # ax11.get_yaxis().set_minor_locator(matplotlib.ticker.AutoMinorLocator())
    ax1.set_xticks([0,25,50,75,100])
    ax1.set_xlim([-3,120])
    ax1.set_xlabel(r"Distance to cell $r$ ($\mu$m)",fontsize=label_fontsize)

    if not normBkgnd:
        ax1.set_ylim([60,250])
        if logScale:
            ax1.set_yticks([60,1e2,1.5e2,2.0e2])
            ax1.set_yticklabels([])
            ax1.set_yticklabels(["0.6","1.0","1.5","2.0"])
            #y_major = matplotlib.ticker.LogLocator(base = 10.0, numticks = 3)
            #ax1.yaxis.set_major_locator(y_major)
            y_minor = matplotlib.ticker.LogLocator(base = 10.0, subs = np.arange(1.0, 10.0) * 0.1, numticks = 10)
            ax1.yaxis.set_minor_locator(y_minor)
            ax1.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
            #ax11.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
        else:
            ax1.set_yticks([50,100,150,200])
            ax1.set_yticklabels([])
            ax1.set_yticklabels(["50","100","150","200"])
    if normBkgnd:
        ax1.set_ylim([0.5,2.7])
        if plotR:
            ax11.set_ylim([0.01,8]) 
   
    
    for i in range(np.shape(NormalOCE_ROI)[0]):
        tmp_df = df_res[df_res["group"]==i]
        ax1.errorbar(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rGmean"],yerr=tmp_df.loc[tmp_df.index[0],"rGcon"],xerr=tmp_df.loc[tmp_df.index[0],"rrcon"],ecolor=colors[1],color=colors[1],capsize=2,zorder=999,linewidth=lw,elinewidth=lw)
        ax1.scatter(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rGmean"],facecolors=colors[1],edgecolors=colors[1],linewidths=lw,s=scatterSize,zorder=1000,marker="o") 
        if plotR:
            ax11.errorbar(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rRmean"],yerr=tmp_df.loc[tmp_df.index[0],"rRcon"],xerr=tmp_df.loc[tmp_df.index[0],"rrcon"],ecolor=colors[1],color=colors[1],capsize=2,linewidth=lw,elinewidth=lw,zorder=80)
            ax11.scatter(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rRmean"],facecolors="white",edgecolors=colors[1],linewidths=lw,s=scatterSize,marker="o",zorder=100) 




    #curve fitting to normal OCE results 
    if curveFit: 
        tmp_df = df_res[df_res["group"]>-1] 
        tmp_df = tmp_df.sort_values("group")
        cx = tmp_df["rrmean"].unique() 
        cy = tmp_df["rGmean"].unique() 
        cx = cx[0:5]
        cy = cy[0:5]
        cy = np.log10(cy)        
        fitx = np.linspace(0.1,70,num=50)
        X = np.log10(cx)
        X = sm.add_constant(X) 
        model = sm.RLM(cy,X) 
        res_fit = model.fit()
        alpha_fit = 0.05
        print(res_fit.summary(alpha=alpha_fit)) 
        print(res_fit.params)
        print(res_fit.conf_int(alpha=alpha_fit))
        p = res_fit.params
        conf = res_fit.conf_int(alpha=alpha_fit)
        fity = p[1] * np.log10(fitx) + p[0] 
        resy = 10**fity 
        #ax1.plot(fitx,resy,color="gray",alpha=0.75,linewidth=2,zorder=0,linestyle="--")
        axCF.annotate(text=r"$\mathrm{G^{\prime}\propto r^{-0.19}}$",xy=(fitx[2],resy[2]),xytext=(5,2.2),arrowprops=dict(facecolor='black',edgecolor="none",width=1.0,headwidth=2,headlength=4,alpha=0.5),zorder=10000,fontsize=8)
        axCF.plot(fitx,resy,color="gray",alpha=0.75,lw=1,zorder=0,linestyle="--")

        for i in range(5):
            tmp_df = df_res[df_res["group"]==i]
            axCF.errorbar(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rGmean"],yerr=tmp_df.loc[tmp_df.index[0],"rGcon"],xerr=tmp_df.loc[tmp_df.index[0],"rrcon"],ecolor=colors[1],color=colors[1],capsize=1,zorder=990,linewidth=lw/1.3,elinewidth=lw/1.3)
            axCF.scatter(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rGmean"],facecolors=colors[1],edgecolors=colors[1],linewidths=lw/1.3,s=scatterSize/5,zorder=995,marker="s") 
       # yci_down = conf[1,0] * np.log10(fitx) + conf[0,0]
       # yci_up = conf[1,1] * np.log10(fitx) + conf[0,1]
       # axCF.fill_between(fitx,10**yci_down,10**yci_up,alpha=0.1,color=colors[1])


    # now including psotCytoD results 
    postCytoD_beadsFileNames = ["pCell19/pBeadsRes19_210405","pCell20/pBeadsRes_210405_6hrs","pCell18/pBeadsRes_210405","pCell17/pBeadsRes_210403","earlyCell11/pCell14/pBeadsRes_210128","earlyCell9/pCell12/pBeadsRes_210124","earlyCell10/pCell13/pBeadsRes_210126"]
    postCytoD_cellFileNames = ["pCell19/pCellBody19_210405","pCell20/pCellBody_210405_6hrs","pCell18/pCellBody_210405","pCell17/pCellBody_210403","earlyCell11/pCell14/pCellBody_210128","earlyCell9/pCell12/pCellBody_210124","earlyCell10/pCell13/pCellBody_210126"]
    postCyotD_N_cell = len(postCytoD_beadsFileNames) 
    postCytoD_xoffset = np.asarray([52,52,52,52,52,52,52,52,52,52,52,52,52,52,52])
    
    postOCE_ROI = np.asarray([[1.2,5.0],[6.1,10],[14,17],[18.1,21],[25,27],[30,35],[40,45],[50,55],[60,65],[70,85],[90,105],[106,115]])
   # postOCE_ROI = np.asarray([[0.0,3.0],[6.1,10],[11.1,14],[18.1,21],[25,27],[30,35],[40,45],[50,55],[60,65],[70,85],[90,105],[106,115]])
    df_res_postCytoD = LoadingFromFileList(root_path,postCytoD_beadsFileNames,postCytoD_cellFileNames,postCytoD_xoffset,CenterToCell=True,normBkgnd=normBkgnd) 
    df_res_postCytoD = SelectROI(postOCE_ROI,df_res_postCytoD,perRange=[10,95])
    df_res_postCytoD.to_csv(os.path.join(save_path,"results_postCytoD.csv"),index=False) 
    print("postCytoD Cell OCE results...")
    print(df_res_postCytoD)  

    
    #df_res_postCytoD["Gre"] = df_res_postCytoD["Gre"] - 20
    for i in range(np.shape(postOCE_ROI)[0]):
        tmp_df = df_res_postCytoD[df_res_postCytoD["group"]==i]
        ax1.errorbar(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rGmean"],yerr=tmp_df.loc[tmp_df.index[0],"rGcon"],xerr=tmp_df.loc[tmp_df.index[0],"rrcon"],ecolor=colors[0],color=colors[0],capsize=2,zorder=990,linewidth=lw,elinewidth=lw)
        ax1.scatter(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rGmean"],facecolors=colors[0],edgecolors=colors[0],linewidths=lw,s=scatterSize,zorder=995,marker="s") 
        if plotR:
            ax11.errorbar(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rRmean"],yerr=tmp_df.loc[tmp_df.index[0],"rRcon"],xerr=tmp_df.loc[tmp_df.index[0],"rrcon"],ecolor=colors[0],color=colors[0],capsize=2,linewidth=lw,elinewidth=lw,zorder=70)
            ax11.scatter(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rRmean"],facecolors="white",edgecolors=colors[0],linewidths=lw,s=scatterSize,marker="s",zorder=90) 

    #curve fitting to postCytoD OCE results 
    if curveFit: 
        tmp_df = df_res_postCytoD[df_res_postCytoD["group"]>-1] 
        tmp_df = tmp_df.sort_values("group")
        cx = tmp_df["rrmean"].unique() 
        cy = tmp_df["rGmean"].unique() 
        cx = cx[0:5]
        cy = cy[0:5]
        cy = np.log10(cy)
        
        fitx = np.linspace(0.1,70,num=50)
        X = np.log10(cx)
        X = sm.add_constant(X) 
        model = sm.RLM(cy,X) 
        res_fit = model.fit()
        print(res_fit.summary(alpha=alpha_fit)) 
        print(res_fit.params)
        print(res_fit.conf_int(alpha=alpha_fit))
        p = res_fit.params
        conf = res_fit.conf_int(alpha=alpha_fit)
        fity = p[1] * np.log10(fitx) + p[0] 
        resy = 10**fity 
        axCF.plot(fitx,resy,color="gray",alpha=0.75,lw=1,zorder=0,linestyle="-.")
        #ax1.plot(fitx,resy,color="gray",alpha=0.75,linewidth=2,zorder=0,linestyle="--")
        axCF.annotate(text=r"$\mathrm{G^{\prime}\propto r^{-0.12}}$",xy=(fitx[1],resy[1]),xytext=(0.7,0.95),arrowprops=dict(facecolor='black',edgecolor="none",width=1.0,headwidth=2,headlength=4,alpha=0.5),zorder=10000,fontsize=8)
        for i in range(5):
            tmp_df = df_res_postCytoD[df_res_postCytoD["group"]==i]
            axCF.errorbar(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rGmean"],yerr=tmp_df.loc[tmp_df.index[0],"rGcon"],xerr=tmp_df.loc[tmp_df.index[0],"rrcon"],ecolor=colors[0],color=colors[0],capsize=1,zorder=990,linewidth=lw/1.3,elinewidth=lw/1.3)
            axCF.scatter(tmp_df.loc[tmp_df.index[0],"rrmean"],tmp_df.loc[tmp_df.index[0],"rGmean"],facecolors=colors[0],edgecolors=colors[0],linewidths=lw/1.3,s=scatterSize/5,zorder=995,marker="s") 
       
      #  yci_down = conf[1,0] * np.log10(fitx) + conf[0,0]
      #  yci_up = conf[1,1] * np.log10(fitx) + conf[0,1]
      #  axCF.fill_between(fitx,10**yci_down,10**yci_up,alpha=0.1,color=colors[0])
    # add legend 
    #fig1.patches.extend([plt.Circle((20,210),4,edgecolor=colors[1],facecolor=colors[1],lw=1,figure=fig1)])   
    #fig1.patches.extend([plt.Circle((20,205),4,edgecolor=colors[1],facecolor="white",lw=1,figure=fig1)])  
    #fig1.patches.extend([plt.Circle((20,200),4,edgecolor=colors[0],facecolor=colors[0],lw=1,figure=fig1)])  
    #fig1.patches.extend([plt.Circle((20,195),4,edgecolor=colors[0],facecolor="white",lw=1,figure=fig1)])  
    ax1.scatter(10,230,s=6,edgecolor=colors[1],facecolor=colors[1],lw=1)
    ax1.scatter(10,220,s=6,edgecolor=colors[1],facecolor="white",lw=1)
    ax1.scatter(10,210,s=6,edgecolor=colors[0],facecolor=colors[0],lw=1)
    ax1.scatter(10,200,s=6,edgecolor=colors[0],facecolor="white",lw=1) 
    
    plt.show()


